#include <iostream>
using namespace std;

struct vec2 {
  float x, y;
};

__global__ void concentracion(float *c, vec2 *linec, float *s, float *tao, int line_size, float deltaL0)
{
  int bx = blockIdx.x;	
  int by = blockIdx.y;
  int tx = threadIdx.x;	
  int ty = threadIdx.y;	
  
  int i = bx*blockDim.x + tx;
  int j = by*blockDim.y + ty;
  
  int Nxf = 512;
  int Nyf = 512;

  float dxf = 1.0 / (Nxf-1);
  float dyf = 1.0 / (Nyf-1);

  vec2 x_im1_n = linec[0]; 

  c[Nxf * j + i] = 0;

/* */
  for (int k=1; k < line_size; k++) {	    	  
    vec2 x_i_n = linec[k];
    vec2 sigma, coord;
    
    sigma.x = x_i_n.x - x_im1_n.x; 
    sigma.y = x_i_n.y - x_im1_n.y;
    float normSig = sqrt ( sigma.x * sigma.x + sigma.y * sigma.y);
    sigma.x /= normSig;
    sigma.y /= normSig;
    coord.x = dxf * i;
    coord.y = dyf * j;
    float midx = 0.5 * ( x_i_n.x + x_im1_n.x ); 
    float midy = 0.5 * ( x_i_n.y + x_im1_n.y );	  
    float aux1 = (coord.x - midx) * sigma.x + (coord.y - midy) * sigma.y;
    float aux2 = (coord.x - midx) * (-sigma.y) + (coord.y - midy) * sigma.x;
    aux1 = aux1 * aux1 / (deltaL0 * deltaL0);
    aux2 = aux2 * aux2 / ( s[k]*s[k] * (1.0 + 4.0 * tao[k]) );
    c[Nxf * j + i] += exp(-aux1-aux2) / (1.7726 * sqrt( 1.0 + 4.0 * tao[k]));
    x_im1_n = linec[k];
  }
/*  */
}	

// Variables _h ya están inicializadas en host
int calcula_concentracion(float *c_h, vec2 *linec_h, float *s_h, float *tao_h, int line_size, float deltaL0)
{
  int size = line_size*sizeof(float);
  vec2 *linec;
  float *s, *tao, *c;

  int Nxf = 512;
  int Nyf = 512;	
  int sizeC = Nxf * Nyf * sizeof(float);

  cudaMalloc((void**) &linec, 2*size);
  cudaMalloc((void**) &s, size);
  cudaMalloc((void**) &tao, size);
  cudaMalloc((void**) &c, sizeC);

/* *
int contadorsote = 0;
cout << endl;
for(int i = 1; i < line_size; i++) {
cout << contadorsote++ << " " << linec_h[i].x << " " << linec_h[i].y << " " << s_h[i] << " " << tao_h[i] << endl;
}
cout << endl; 
int pausota; cin >> pausota;
* */

  cudaMemcpy(linec, linec_h, 2*size, cudaMemcpyHostToDevice);
  cudaMemcpy(s, s_h, size, cudaMemcpyHostToDevice);
  cudaMemcpy(tao, tao_h, size, cudaMemcpyHostToDevice);
  dim3 DimGrid(32, 32, 1);
  dim3 DimBlock(16, 16, 1);
  
  concentracion<<<DimGrid, DimBlock>>>(c, linec, s, tao, line_size, deltaL0);
  
  cudaMemcpy(c_h, c, sizeC, cudaMemcpyDeviceToHost ); 

/* *
      cout << endl;

        for (int j = 0; j < Nyf; j++) {
      for(int i  = 0; i < Nxf; i++) {
	  cout << c_h[i*Nxf+j] << " ";
	}
	cout << endl;
      }
      int pausota; cin >> pausota;
/* */

  cudaFree(linec);
  cudaFree(s);
  cudaFree(tao);
  cudaFree(c);

  return 0;
}

